Differential abundance analysis

In [1]:
rm(list=ls())
# Set the working directory to where folders named after the samples are located. 
# The folder contains spliced.mtx, unspliced.mtx, barcodes and gene id files, and json files produced by scKB that documents the sequencing stats. 
setwd("/scratch/AG_Ohler/CheWei/scKB")
In [1]:
# Load libraries
library(edgeR)
library(tidyverse)
library(Seurat)
library(cowplot)
library(ComplexHeatmap)
library(circlize)
library(GeneOverlap)
library(gprofiler2)
library(ggrepel)
library(ggplot2)
library(scater)
library(RColorBrewer)
Loading required package: limma

── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──

 ggplot2 3.3.0      purrr   0.3.3
 tibble  3.0.1      dplyr   0.8.5
 tidyr   1.0.2      stringr 1.4.0
 readr   1.3.1      forcats 0.5.0

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
 dplyr::filter() masks stats::filter()
 dplyr::lag()    masks stats::lag()


********************************************************

Note: As of version 1.0.0, cowplot does not change the

  default ggplot2 theme anymore. To recover the previous

  behavior, execute:
  theme_set(theme_cowplot())

********************************************************


Loading required package: grid

========================================
ComplexHeatmap version 2.2.0
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
  genomic data. Bioinformatics 2016.
========================================


========================================
circlize version 0.4.8
CRAN page: https://cran.r-project.org/package=circlize
Github page: https://github.com/jokergoo/circlize
Documentation: http://jokergoo.github.io/circlize_book/book/

If you use it in published research, please cite:
Gu, Z. circlize implements and enhances circular visualization 
  in R. Bioinformatics 2014.
========================================


Loading required package: SingleCellExperiment

Loading required package: SummarizedExperiment

Loading required package: GenomicRanges

Loading required package: stats4

Loading required package: BiocGenerics

Loading required package: parallel


Attaching package: ‘BiocGenerics’


The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB


The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union


The following object is masked from ‘package:limma’:

    plotMA


The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs


The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames,
    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which, which.max, which.min


Loading required package: S4Vectors


Attaching package: ‘S4Vectors’


The following objects are masked from ‘package:dplyr’:

    first, rename


The following object is masked from ‘package:tidyr’:

    expand


The following object is masked from ‘package:base’:

    expand.grid


Loading required package: IRanges


Attaching package: ‘IRanges’


The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice


The following object is masked from ‘package:purrr’:

    reduce


Loading required package: GenomeInfoDb

Loading required package: Biobase

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.


Loading required package: DelayedArray

Loading required package: matrixStats


Attaching package: ‘matrixStats’


The following objects are masked from ‘package:Biobase’:

    anyMissing, rowMedians


The following object is masked from ‘package:dplyr’:

    count


Loading required package: BiocParallel


Attaching package: ‘DelayedArray’


The following objects are masked from ‘package:matrixStats’:

    colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges


The following object is masked from ‘package:purrr’:

    simplify


The following objects are masked from ‘package:base’:

    aperm, apply, rowsum



Attaching package: ‘SummarizedExperiment’


The following object is masked from ‘package:Seurat’:

    Assays



Attaching package: ‘SingleCellExperiment’


The following object is masked from ‘package:edgeR’:

    cpm



Attaching package: ‘scater’


The following object is masked from ‘package:limma’:

    plotMDS


In [2]:
library(future)
#for 300gb ram 
options(future.globals.maxSize = 300000 * 1024^2)
Attaching package: ‘future’


The following object is masked from ‘package:SummarizedExperiment’:

    values


The following object is masked from ‘package:GenomicRanges’:

    values


The following object is masked from ‘package:IRanges’:

    values


The following object is masked from ‘package:S4Vectors’:

    values


In [3]:
rc.integrated <- readRDS("./rc.integrated_wt_shr_scr.rds")
In [4]:
feature_names <- read_tsv("./supp_data/features.tsv.gz", col_names = c("AGI", "Name", "Type")) %>%
  select(-Type) %>%
  distinct()
Parsed with column specification:
cols(
  AGI = col_character(),
  Name = col_character(),
  Type = col_character()
)

1. Add metadata and organize by genotype and experiment date

In [5]:
bscs <- read.csv("./Benfey_single_cell-Samples.csv", na.strings=c("","NA"), stringsAsFactors = F)
  bscs <- bscs %>% select(c('sample','name','source','geno', 'genotype','transgene','treatment','age','timepoint','rep','target_cells','date','seq_run'))
  bscs$date <- gsub('^([0-9]{4})([0-9]{2})([0-9]+)$', '\\1-\\2-\\3', bscs$date)
  bscs$target_cells <- prettyNum(bscs$target_cells, big.mark = ',')

orig.idents <- data.frame(sample=rc.integrated$orig.ident)

bscs %>% filter(sample %in% orig.idents$sample)
A data.frame: 10 × 13
samplenamesourcegenogenotypetransgenetreatmentagetimepointreptarget_cellsdateseq_run
<chr><chr><chr><chr><chr><chr><chr><chr><dbl><int><chr><chr><chr>
sc_12WT Col-0 untreatedBenfey labWT WT Col-0 NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_20WT Col-0_RS1 Benfey labWT WT Col-0 NAUntreated5_dayNANA10,000NA Shahan_6177
sc_21WT Col-0_RS2 Benfey labWT WT Col-0 NAUntreated5_dayNANA10,000NA Shahan_6177
sc_25scr-4_1 Benfey labscr-4scr-4 (Ws backcrossed to Col?)NAUntreated5_dayNANA10,000NA Shahan_6177
sc_30WT Col-0_RS3 Benfey labWT WT Col-0 NAUntreated5_dayNANA10,0002020-02-10Nolan_6199 (NextSeq); Nolan_6226 (NovaSeq S4)
sc_31WT Col-0_RS4 Benfey labWT WT Col-0 NAUntreated5_dayNANA10,0002020-02-10Nolan_6199 (NextSeq); Nolan_6226 (NovaSeq S4)
sc_36scr_4_2 Benfey labscr-4scr-4 NAUntreated5_dayNANA10,0002020-02-10Nolan_6199 (NextSeq); Nolan_6226 (NovaSeq S4)
sc_51WT Col (RS_5) Benfey labWT WT Col-0 NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_52shr-2 Benfey labshr-2shr-2 NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2 Benfey labshr-2shr-2 NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
In [6]:
orig.idents_info <- left_join(orig.idents, bscs, by="sample")
orig.idents_info
Warning message:
“Column `sample` joining factor and character vector, coercing into character vector”
A data.frame: 86260 × 13
samplenamesourcegenogenotypetransgenetreatmentagetimepointreptarget_cellsdateseq_run
<chr><chr><chr><chr><chr><chr><chr><chr><dbl><int><chr><chr><chr>
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_12WT Col-0 untreatedBenfey labWTWT Col-0NAUntreated5_dayNANA10,0002019-12-20Nolan_6131;Shahan_6158
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
sc_53shr-2Benfey labshr-2shr-2NAUntreated5_dayNANA10,0002020-02-24Nolan_6226
In [7]:
treatment <- orig.idents_info$treatment
names(treatment) <- colnames(rc.integrated)
head(treatment)
rc.integrated <- AddMetaData(
  object = rc.integrated,
  metadata = treatment,
  col.name = 'treatment')


timepoint <- orig.idents_info$timepoint
names(timepoint) <- colnames(rc.integrated)

rc.integrated <- AddMetaData(
  object = rc.integrated,
  metadata = timepoint,
  col.name = 'timepoint')



age <- orig.idents_info$age
names(age) <- colnames(rc.integrated)
rc.integrated <- AddMetaData(
  object = rc.integrated,
  metadata = age,
  col.name = 'age')



rep <- orig.idents_info$rep
names(rep) <- colnames(rc.integrated)

rc.integrated <- AddMetaData(
  object = rc.integrated,
  metadata = rep,
  col.name = 'rep')

date <- orig.idents_info$date
names(date) <- colnames(rc.integrated)

rc.integrated <- AddMetaData(
  object = rc.integrated,
  metadata = date,
  col.name = 'date')

geno <- orig.idents_info$geno

names(geno) <- colnames(rc.integrated)

rc.integrated <- AddMetaData(
  object = rc.integrated,
  metadata = geno,
  col.name = 'geno')
AAACCCAAGCCATTCA_1
'Untreated'
AAACCCAAGCTGTTCA_1
'Untreated'
AAACCCAAGGCGAACT_1
'Untreated'
AAACCCACACCAGCCA_1
'Untreated'
AAACCCACACTGGCGT_1
'Untreated'
AAACCCACAGACCTGC_1
'Untreated'
In [8]:
table(rc.integrated$treatment)
Untreated 
    86260 
In [9]:
table(rc.integrated$geno)
scr-4 shr-2    WT 
13742 14612 57906 
In [10]:
rc.integrated$geno_trt <- rep("WT")
rc.integrated$geno_trt[rc.integrated$geno=="WT"] <- "WT"
rc.integrated$geno_trt[rc.integrated$geno=="shr-2"] <- "shr"
rc.integrated$geno_trt[rc.integrated$geno=="scr-4"] <- "scr"
rc.integrated$geno_trt <- factor(rc.integrated$geno_trt, levels=c("WT", "shr", "scr"))
In [11]:
table(rc.integrated$geno_trt)
   WT   shr   scr 
57906 14612 13742 
In [12]:
table(rc.integrated$celltype.anno)
Putative Quiescent Center           Stem Cell Niche                 Columella 
                      397                       831                     10548 
         Lateral Root Cap              Atrichoblast               Trichoblast 
                    13856                     11362                     10742 
                   Cortex                Endodermis                 Pericycle 
                     9719                      6064                      6379 
                   Phloem                     Xylem                Procambium 
                     4095                      3964                      8303 
In [13]:
# Simple QC label
rc.integrated$celltype.anno <- gsub("Putative Quiescent Center", "Quiescent Center", rc.integrated$celltype.anno, ignore.case = FALSE, perl = FALSE,
     fixed = T, useBytes = FALSE)

order <- c("Quiescent Center", "Stem Cell Niche", "Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Pericycle", "Phloem", "Xylem", "Procambium", "Unknown")
palette <- c("#9400D3","#DCD0FF", "#5AB953", "#BFEF45", "#008080", "#21B6A8", "#82B6FF", "#0000FF","#FF9900","#E6194B", "#9A6324", "#FFE119","#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno, levels = order[sort(match(unique(rc.integrated$celltype.anno),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
In [14]:
options(repr.plot.width=6, repr.plot.height=6)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols=color)
In [15]:
options(repr.plot.width=12, repr.plot.height=6)
DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", split.by="geno_trt", cols=color)

2. Draw a heatmap of each celltype's numbers

In [16]:
sample_num_df <- data.frame(table(rc.integrated$geno_trt))

sample_num_df <- dplyr::rename(sample_num_df, "total_cells"="Freq", "orig.id"="Var1")

sample_num_df

rc.integrated$sample.cell.p <- paste(rc.integrated$geno_trt, rc.integrated$celltype.anno, sep = ",")

sample_cell_df <- data.frame(table(rc.integrated$sample.cell.p))

sample_cell_df <- separate(sample_cell_df, col="Var1", sep=",", into=c("orig.id", "celltype")) %>%
  dplyr::rename("n_celltype"="Freq") %>%
  left_join(sample_num_df) %>%
  mutate(proportion=n_celltype/total_cells)

sample_cell_df

wt_cells <- filter(sample_cell_df, orig.id == "WT") %>%
  group_by(celltype) %>%
  summarise(wt_prop=mean(proportion))

wt_cells

sample_cell_df <- left_join(sample_cell_df, wt_cells) %>%
  mutate(wt_norm_prop=proportion/wt_prop)

sample_cell_df


sample_cell_wide <- sample_cell_df %>%
  ungroup() %>%
  select(orig.id, celltype, wt_norm_prop) %>%
  pivot_wider(names_from = celltype, values_from = wt_norm_prop)

sample_cell_wide

sample_cell_wide[is.na(sample_cell_wide)==T] <- 0.0000001

sample_cell_wide

cell_m <- as.matrix(log2(sample_cell_wide[, -1]))
rownames(cell_m) <- sample_cell_wide$orig.id

cell_m
A data.frame: 3 × 2
orig.idtotal_cells
<fct><int>
WT 57906
shr14612
scr13742
Joining, by = "orig.id"

Warning message:
“Column `orig.id` joining character vector and factor, coercing into character vector”
A data.frame: 36 × 5
orig.idcelltypen_celltypetotal_cellsproportion
<chr><chr><int><int><dbl>
scrAtrichoblast 2009137420.146194149
scrColumella 2161137420.157255130
scrCortex 946137420.068840052
scrEndodermis 361137420.026269830
scrLateral Root Cap3637137420.264663077
scrPericycle 453137420.032964634
scrPhloem 566137420.041187600
scrProcambium 986137420.071750837
scrQuiescent Center 81137420.005894339
scrStem Cell Niche 166137420.012079755
scrTrichoblast 1951137420.141973512
scrXylem 425137420.030927085
shrAtrichoblast 2064146120.141253764
shrColumella 2888146120.197645771
shrCortex 2265146120.155009581
shrEndodermis 190146120.013003011
shrLateral Root Cap3249146120.222351492
shrPericycle 318146120.021762935
shrPhloem 132146120.009033671
shrProcambium 1463146120.100123186
shrQuiescent Center 73146120.004995894
shrStem Cell Niche 268146120.018341090
shrTrichoblast 1399146120.095743225
shrXylem 303146120.020736381
WT Atrichoblast 7289579060.125876420
WT Columella 5499579060.094964252
WT Cortex 6508579060.112389044
WT Endodermis 5513579060.095206024
WT Lateral Root Cap6970579060.120367492
WT Pericycle 5608579060.096846613
WT Phloem 3397579060.058664042
WT Procambium 5854579060.101094878
WT Quiescent Center 243579060.004196456
WT Stem Cell Niche 397579060.006855939
WT Trichoblast 7392579060.127655165
WT Xylem 3236579060.055883674
A tibble: 12 × 2
celltypewt_prop
<chr><dbl>
Atrichoblast 0.125876420
Columella 0.094964252
Cortex 0.112389044
Endodermis 0.095206024
Lateral Root Cap0.120367492
Pericycle 0.096846613
Phloem 0.058664042
Procambium 0.101094878
Quiescent Center0.004196456
Stem Cell Niche 0.006855939
Trichoblast 0.127655165
Xylem 0.055883674
Joining, by = "celltype"

A data.frame: 36 × 7
orig.idcelltypen_celltypetotal_cellsproportionwt_propwt_norm_prop
<chr><chr><int><int><dbl><dbl><dbl>
scrAtrichoblast 2009137420.1461941490.1258764201.1614101
scrColumella 2161137420.1572551300.0949642521.6559403
scrCortex 946137420.0688400520.1123890440.6125157
scrEndodermis 361137420.0262698300.0952060240.2759261
scrLateral Root Cap3637137420.2646630770.1203674922.1987920
scrPericycle 453137420.0329646340.0968466130.3403798
scrPhloem 566137420.0411876000.0586640420.7020928
scrProcambium 986137420.0717508370.1010948780.7097376
scrQuiescent Center 81137420.0058943390.0041964561.4045990
scrStem Cell Niche 166137420.0120797550.0068559391.7619404
scrTrichoblast 1951137420.1419735120.1276551651.1121643
scrXylem 425137420.0309270850.0558836740.5534190
shrAtrichoblast 2064146120.1412537640.1258764201.1221622
shrColumella 2888146120.1976457710.0949642522.0812650
shrCortex 2265146120.1550095810.1123890441.3792232
shrEndodermis 190146120.0130030110.0952060240.1365776
shrLateral Root Cap3249146120.2223514920.1203674921.8472719
shrPericycle 318146120.0217629350.0968466130.2247155
shrPhloem 132146120.0090336710.0586640420.1539899
shrProcambium 1463146120.1001231860.1010948780.9903883
shrQuiescent Center 73146120.0049958940.0041964561.1905030
shrStem Cell Niche 268146120.0183410900.0068559392.6752119
shrTrichoblast 1399146120.0957432250.1276551650.7500145
shrXylem 303146120.0207363810.0558836740.3710633
WT Atrichoblast 7289579060.1258764200.1258764201.0000000
WT Columella 5499579060.0949642520.0949642521.0000000
WT Cortex 6508579060.1123890440.1123890441.0000000
WT Endodermis 5513579060.0952060240.0952060241.0000000
WT Lateral Root Cap6970579060.1203674920.1203674921.0000000
WT Pericycle 5608579060.0968466130.0968466131.0000000
WT Phloem 3397579060.0586640420.0586640421.0000000
WT Procambium 5854579060.1010948780.1010948781.0000000
WT Quiescent Center 243579060.0041964560.0041964561.0000000
WT Stem Cell Niche 397579060.0068559390.0068559391.0000000
WT Trichoblast 7392579060.1276551650.1276551651.0000000
WT Xylem 3236579060.0558836740.0558836741.0000000
A tibble: 3 × 13
orig.idAtrichoblastColumellaCortexEndodermisLateral Root CapPericyclePhloemProcambiumQuiescent CenterStem Cell NicheTrichoblastXylem
<chr><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
scr1.1614101.6559400.61251570.27592612.1987920.34037980.70209280.70973761.4045991.7619401.11216430.5534190
shr1.1221622.0812651.37922320.13657761.8472720.22471550.15398990.99038831.1905032.6752120.75001450.3710633
WT 1.0000001.0000001.00000001.00000001.0000001.00000001.00000001.00000001.0000001.0000001.00000001.0000000
A tibble: 3 × 13
orig.idAtrichoblastColumellaCortexEndodermisLateral Root CapPericyclePhloemProcambiumQuiescent CenterStem Cell NicheTrichoblastXylem
<chr><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
scr1.1614101.6559400.61251570.27592612.1987920.34037980.70209280.70973761.4045991.7619401.11216430.5534190
shr1.1221622.0812651.37922320.13657761.8472720.22471550.15398990.99038831.1905032.6752120.75001450.3710633
WT 1.0000001.0000001.00000001.00000001.0000001.00000001.00000001.00000001.0000001.0000001.00000001.0000000
A matrix: 3 × 12 of type dbl
AtrichoblastColumellaCortexEndodermisLateral Root CapPericyclePhloemProcambiumQuiescent CenterStem Cell NicheTrichoblastXylem
scr0.21587750.7276506-0.7071813-1.8576461.1367111-1.554783-0.5102664-0.494642340.49015840.8171651 0.1533699-0.853556
shr0.16628131.0574606 0.4638560-2.8722070.8853963-2.153828-2.6990922-0.013933790.25157121.4196532-0.4150096-1.430263
WT0.00000000.0000000 0.0000000 0.0000000.0000000 0.000000 0.0000000 0.000000000.00000000.0000000 0.0000000 0.000000
In [17]:
# add to hm cell_fun = function(j, i, x, y, width, height, fill) {grid.text(sprintf("%.0f", int_matrix[i, j]), x, y, gp = gpar(fontsize = 10))})

Cell_hm <- Heatmap(cell_m, name = "log2(mutant/WT cell proportion)", 
                   heatmap_legend_param = list(title_position="topcenter", 
                                               color_bar = "continuous", 
                                               legend_direction = "horizontal"), 
                   col = colorRamp2(c(-3, 0, 3), c("blue", "beige", "#e31a1c")), 
                   cluster_rows = F, 
                   cluster_columns = T, 
                   use_raster= FALSE, 
                   show_column_names = TRUE, show_row_names = TRUE, show_row_dend = TRUE, show_column_dend = TRUE,
                   clustering_distance_rows = "pearson", clustering_distance_columns = "pearson",
                 row_names_gp = gpar(fontsize = 12)) 


draw(Cell_hm, padding = unit(c(10, 10, 10, 10), "mm"), heatmap_legend_side = "top")
In [18]:
sample_num_df <- data.frame(table(rc.integrated$orig.ident))

sample_num_df <- dplyr::rename(sample_num_df, "total_cells"="Freq", "orig.id"="Var1")

sample_num_df

rc.integrated$sample.cell.p <- paste(rc.integrated$orig.ident, rc.integrated$celltype.anno, sep = ",")

sample_cell_df <- data.frame(table(rc.integrated$sample.cell.p))

sample_cell_df <- separate(sample_cell_df, col="Var1", sep=",", into=c("orig.id", "celltype")) %>%
  dplyr::rename("n_celltype"="Freq") %>%
  left_join(sample_num_df) %>%
  mutate(proportion=n_celltype/total_cells)

sample_cell_df

wt_cells <- filter(sample_cell_df, orig.id %in% c("sc_12", "sc_20", "sc_21", "sc_30", "sc_31", "sc_51")) %>%
  group_by(celltype) %>%
  summarise(wt_prop=mean(proportion))

wt_cells

sample_cell_df <- left_join(sample_cell_df, wt_cells) %>%
  mutate(wt_norm_prop=proportion/wt_prop)

sample_cell_df


sample_cell_wide <- sample_cell_df %>%
  ungroup() %>%
  select(orig.id, celltype, wt_norm_prop) %>%
  pivot_wider(names_from = celltype, values_from = wt_norm_prop)

sample_cell_wide

sample_cell_wide[is.na(sample_cell_wide)==T] <- 0.0000001

sample_cell_wide

cell_m <- as.matrix(log2(sample_cell_wide[, -1]))
rownames(cell_m) <- sample_cell_wide$orig.id

cell_m
A data.frame: 10 × 2
orig.idtotal_cells
<fct><int>
sc_1210293
sc_2012591
sc_21 9658
sc_25 7736
sc_3010195
sc_31 8775
sc_36 6006
sc_51 6394
sc_52 6241
sc_53 8371
Joining, by = "orig.id"

Warning message:
“Column `orig.id` joining character vector and factor, coercing into character vector”
A data.frame: 120 × 5
orig.idcelltypen_celltypetotal_cellsproportion
<chr><chr><int><int><dbl>
sc_12Atrichoblast 1351102930.131254250
sc_12Columella 1309102930.127173807
sc_12Cortex 1146102930.111337802
sc_12Endodermis 664102930.064509861
sc_12Lateral Root Cap1328102930.129019722
sc_12Pericycle 1159102930.112600797
sc_12Phloem 410102930.039832896
sc_12Procambium 1096102930.106480132
sc_12Quiescent Center 48102930.004663363
sc_12Stem Cell Niche 65102930.006314971
sc_12Trichoblast 1313102930.127562421
sc_12Xylem 404102930.039249976
sc_20Atrichoblast 1556125910.123580335
sc_20Columella 1116125910.088634739
sc_20Cortex 1404125910.111508220
sc_20Endodermis 947125910.075212453
sc_20Lateral Root Cap1446125910.114843936
sc_20Pericycle 1153125910.091573346
sc_20Phloem 818125910.064967040
sc_20Procambium 1324125910.105154475
sc_20Quiescent Center 50125910.003971090
sc_20Stem Cell Niche 91125910.007227385
sc_20Trichoblast 1865125910.148121674
sc_20Xylem 821125910.065205305
sc_21Atrichoblast 1120 96580.115966039
sc_21Columella 608 96580.062952992
sc_21Cortex 1081 96580.111927935
sc_21Endodermis 916 96580.094843653
sc_21Lateral Root Cap 935 96580.096810934
sc_21Pericycle 1100 96580.113895216
sc_51Phloem 33663940.052549265
sc_51Procambium 55063940.086018142
sc_51Quiescent Center 5563940.008601814
sc_51Stem Cell Niche 7063940.010947764
sc_51Trichoblast 63363940.098999062
sc_51Xylem 28863940.045042227
sc_52Atrichoblast 76262410.122095818
sc_52Columella 132262410.211825028
sc_52Cortex 100162410.160390963
sc_52Endodermis 9562410.015221920
sc_52Lateral Root Cap138362410.221599103
sc_52Pericycle 13562410.021631149
sc_52Phloem 4662410.007370614
sc_52Procambium 63062410.100945361
sc_52Quiescent Center 3262410.005127383
sc_52Stem Cell Niche 10162410.016183304
sc_52Trichoblast 61962410.099182823
sc_52Xylem 11562410.018426534
sc_53Atrichoblast 130283710.155536973
sc_53Columella 156683710.187074424
sc_53Cortex 126483710.150997491
sc_53Endodermis 9583710.011348704
sc_53Lateral Root Cap186683710.222912436
sc_53Pericycle 18383710.021861187
sc_53Phloem 8683710.010273563
sc_53Procambium 83383710.099510214
sc_53Quiescent Center 4183710.004897862
sc_53Stem Cell Niche 16783710.019949827
sc_53Trichoblast 78083710.093178832
sc_53Xylem 18883710.022458488
A tibble: 12 × 2
celltypewt_prop
<chr><dbl>
Atrichoblast 0.124761278
Columella 0.096459978
Cortex 0.112091770
Endodermis 0.096918446
Lateral Root Cap0.123867450
Pericycle 0.096059388
Phloem 0.058352882
Procambium 0.100091969
Quiescent Center0.004538212
Stem Cell Niche 0.007126989
Trichoblast 0.124699541
Xylem 0.055032097
Joining, by = "celltype"

A data.frame: 120 × 7
orig.idcelltypen_celltypetotal_cellsproportionwt_propwt_norm_prop
<chr><chr><int><int><dbl><dbl><dbl>
sc_12Atrichoblast 1351102930.1312542500.1247612781.0520432
sc_12Columella 1309102930.1271738070.0964599781.3184101
sc_12Cortex 1146102930.1113378020.1120917700.9932737
sc_12Endodermis 664102930.0645098610.0969184460.6656097
sc_12Lateral Root Cap1328102930.1290197220.1238674501.0415950
sc_12Pericycle 1159102930.1126007970.0960593881.1721998
sc_12Phloem 410102930.0398328960.0583528820.6826209
sc_12Procambium 1096102930.1064801320.1000919691.0638229
sc_12Quiescent Center 48102930.0046633630.0045382121.0275773
sc_12Stem Cell Niche 65102930.0063149710.0071269890.8860644
sc_12Trichoblast 1313102930.1275624210.1246995411.0229582
sc_12Xylem 404102930.0392499760.0550320970.7132197
sc_20Atrichoblast 1556125910.1235803350.1247612780.9905344
sc_20Columella 1116125910.0886347390.0964599780.9188758
sc_20Cortex 1404125910.1115082200.1120917700.9947940
sc_20Endodermis 947125910.0752124530.0969184460.7760386
sc_20Lateral Root Cap1446125910.1148439360.1238674500.9271519
sc_20Pericycle 1153125910.0915733460.0960593880.9532993
sc_20Phloem 818125910.0649670400.0583528821.1133476
sc_20Procambium 1324125910.1051544750.1000919691.0505786
sc_20Quiescent Center 50125910.0039710900.0045382120.8750342
sc_20Stem Cell Niche 91125910.0072273850.0071269891.0140867
sc_20Trichoblast 1865125910.1481216740.1246995411.1878285
sc_20Xylem 821125910.0652053050.0550320971.1848595
sc_21Atrichoblast 1120 96580.1159660390.1247612780.9295035
sc_21Columella 608 96580.0629529920.0964599780.6526333
sc_21Cortex 1081 96580.1119279350.1120917700.9985384
sc_21Endodermis 916 96580.0948436530.0969184460.9785924
sc_21Lateral Root Cap 935 96580.0968109340.1238674500.7815688
sc_21Pericycle 1100 96580.1138952160.0960593881.1856750
sc_51Phloem 33663940.0525492650.0583528820.9005428
sc_51Procambium 55063940.0860181420.1000919690.8593911
sc_51Quiescent Center 5563940.0086018140.0045382121.8954192
sc_51Stem Cell Niche 7063940.0109477640.0071269891.5360994
sc_51Trichoblast 63363940.0989990620.1246995410.7939008
sc_51Xylem 28863940.0450422270.0550320970.8184719
sc_52Atrichoblast 76262410.1220958180.1247612780.9786355
sc_52Columella 132262410.2118250280.0964599782.1959888
sc_52Cortex 100162410.1603909630.1120917701.4308897
sc_52Endodermis 9562410.0152219200.0969184460.1570591
sc_52Lateral Root Cap138362410.2215991030.1238674501.7890019
sc_52Pericycle 13562410.0216311490.0960593880.2251852
sc_52Phloem 4662410.0073706140.0583528820.1263110
sc_52Procambium 63062410.1009453610.1000919691.0085261
sc_52Quiescent Center 3262410.0051273830.0045382121.1298246
sc_52Stem Cell Niche 10162410.0161833040.0071269892.2707071
sc_52Trichoblast 61962410.0991828230.1246995410.7953744
sc_52Xylem 11562410.0184265340.0550320970.3348325
sc_53Atrichoblast 130283710.1555369730.1247612781.2466767
sc_53Columella 156683710.1870744240.0964599781.9393994
sc_53Cortex 126483710.1509974910.1120917701.3470881
sc_53Endodermis 9583710.0113487040.0969184460.1170954
sc_53Lateral Root Cap186683710.2229124360.1238674501.7996046
sc_53Pericycle 18383710.0218611870.0960593880.2275799
sc_53Phloem 8683710.0102735630.0583528820.1760592
sc_53Procambium 83383710.0995102140.1000919690.9941878
sc_53Quiescent Center 4183710.0048978620.0045382121.0792492
sc_53Stem Cell Niche 16783710.0199498270.0071269892.7991944
sc_53Trichoblast 78083710.0931788320.1246995410.7472267
sc_53Xylem 18883710.0224584880.0550320970.4080980
A tibble: 10 × 13
orig.idAtrichoblastColumellaCortexEndodermisLateral Root CapPericyclePhloemProcambiumQuiescent CenterStem Cell NicheTrichoblastXylem
<chr><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
sc_121.05204321.31841010.99327370.66560971.04159501.17219980.68262091.06382291.02757730.88606441.02295820.7132197
sc_200.99053440.91887580.99479400.77603860.92715190.95329931.11334761.05057860.87503421.01408671.18782851.1848595
sc_210.92950350.65263330.99853840.97859240.78156881.18567501.04511911.17411170.15970780.61017731.20978301.2568203
sc_251.19773731.64429960.77611290.28142302.22700280.33238440.69337080.64831791.28177351.01570010.92258990.6623943
sc_301.18323080.88060981.00457171.37741390.85680670.91389431.05730700.86041520.25936370.67437700.96121190.9767362
sc_310.96092350.97940031.08376811.12762590.96233740.92179441.20106270.99168061.78289781.27919510.82431761.0498923
sc_361.13837121.61218320.40551190.25769112.02030280.35706070.72189310.80512031.32078582.56981161.41665860.4326479
sc_510.88376471.25007070.92505421.07471961.43054020.85313710.90054280.85939111.89541921.53609940.79390080.8184719
sc_520.97863552.19598881.43088970.15705911.78900190.22518520.12631101.00852611.12982462.27070710.79537440.3348325
sc_531.24667671.93939941.34708810.11709541.79960460.22757990.17605920.99418781.07924922.79919440.74722670.4080980
A tibble: 10 × 13
orig.idAtrichoblastColumellaCortexEndodermisLateral Root CapPericyclePhloemProcambiumQuiescent CenterStem Cell NicheTrichoblastXylem
<chr><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
sc_121.05204321.31841010.99327370.66560971.04159501.17219980.68262091.06382291.02757730.88606441.02295820.7132197
sc_200.99053440.91887580.99479400.77603860.92715190.95329931.11334761.05057860.87503421.01408671.18782851.1848595
sc_210.92950350.65263330.99853840.97859240.78156881.18567501.04511911.17411170.15970780.61017731.20978301.2568203
sc_251.19773731.64429960.77611290.28142302.22700280.33238440.69337080.64831791.28177351.01570010.92258990.6623943
sc_301.18323080.88060981.00457171.37741390.85680670.91389431.05730700.86041520.25936370.67437700.96121190.9767362
sc_310.96092350.97940031.08376811.12762590.96233740.92179441.20106270.99168061.78289781.27919510.82431761.0498923
sc_361.13837121.61218320.40551190.25769112.02030280.35706070.72189310.80512031.32078582.56981161.41665860.4326479
sc_510.88376471.25007070.92505421.07471961.43054020.85313710.90054280.85939111.89541921.53609940.79390080.8184719
sc_520.97863552.19598881.43088970.15705911.78900190.22518520.12631101.00852611.12982462.27070710.79537440.3348325
sc_531.24667671.93939941.34708810.11709541.79960460.22757990.17605920.99418781.07924922.79919440.74722670.4080980
A matrix: 10 × 12 of type dbl
AtrichoblastColumellaCortexEndodermisLateral Root CapPericyclePhloemProcambiumQuiescent CenterStem Cell NicheTrichoblastXylem
sc_12 0.07319391 0.39879918-0.009736839-0.58725157 0.05879449 0.22921851-0.55084354 0.089258051 0.03924689-0.17451647 0.03274722-0.48758155
sc_20-0.01372104-0.12205823-0.007530284-0.36579972-0.10912244-0.06899888 0.15490406 0.071184037-0.19258875 0.02018102 0.24832660 0.24471604
sc_21-0.10546787-0.61565553-0.002110197-0.03122005-0.35555522 0.24570863 0.06366735 0.231569713-2.64649352-0.71269969 0.27474835 0.32977842
sc_25 0.26031147 0.71747316-0.365661641-1.82918798 1.15510335-1.58907528-0.52830096-0.625226647 0.35814132 0.02247447-0.11623862-0.59423785
sc_30 0.24273148-0.18342523 0.006580516 0.46196211-0.22295840-0.12990071 0.08039433-0.216895146-1.94695135-0.56837268-0.05707365-0.03395919
sc_31-0.05750656-0.03002941 0.116056057 0.17328850-0.05538530-0.11748305 0.26431147-0.012052623 0.83422401 0.35523632-0.27872785 0.07024141
sc_36 0.18697104 0.68901569-1.302183875-1.95628516 1.01457150-1.48575866-0.47014291-0.312723643 0.40139655 1.36166259 0.50249212-1.20873454
sc_51-0.17826572 0.32200973-0.112390225 0.10396025 0.51656009-0.22915046-0.15113333-0.218613341 0.92251699 0.61927162-0.33296940-0.28899515
sc_52-0.03115645 1.13487067 0.516912506-2.67062098 0.83915491-2.15081633-2.98494726 0.012248401 0.17609881 1.18314162-0.33029396-1.57848857
sc_53 0.31808733 0.95560994 0.429844230-3.09424374 0.84767998-2.13555485-2.50586723-0.008409698 0.11002805 1.48501166-0.42038201-1.29301253
In [19]:
options(repr.plot.width=12, repr.plot.height=10)
# add to hm cell_fun = function(j, i, x, y, width, height, fill) {grid.text(sprintf("%.0f", int_matrix[i, j]), x, y, gp = gpar(fontsize = 10))})

Cell_hm <- Heatmap(cell_m, name = "log2(sample/WT cell proportion)", 
                   heatmap_legend_param = list(title_position="topcenter", 
                                               color_bar = "continuous", 
                                               legend_direction = "horizontal"), 
                   col = colorRamp2(c(-2, 0, 2), c("blue", "beige", "#e31a1c")), 
                   cluster_rows = T, 
                   cluster_columns = T, 
                   use_raster= FALSE, 
                   show_column_names = TRUE, show_row_names = TRUE, show_row_dend = TRUE, show_column_dend = TRUE,
                   clustering_distance_rows = "pearson", clustering_distance_columns = "pearson",
                 row_names_gp = gpar(fontsize = 12)) 


draw(Cell_hm, padding = unit(c(10, 10, 10, 10), "mm"), heatmap_legend_side = "top")
In [20]:
table(rc.integrated$batch)
batch 1 batch 2 batch 3 batch 4 
  10293   29985   24976   21006 
In [21]:
table(rc.integrated$batch, rc.integrated$orig.ident)
         
          sc_12 sc_20 sc_21 sc_25 sc_30 sc_31 sc_36 sc_51 sc_52 sc_53
  batch 1 10293     0     0     0     0     0     0     0     0     0
  batch 2     0 12591  9658  7736     0     0     0     0     0     0
  batch 3     0     0     0     0 10195  8775  6006     0     0     0
  batch 4     0     0     0     0     0     0     0  6394  6241  8371
In [22]:
rc.integrated$batch<- gsub(" ", "_", rc.integrated$batch, ignore.case = FALSE, perl = FALSE,
     fixed = T, useBytes = FALSE)
table(rc.integrated$batch, rc.integrated$orig.ident)
         
          sc_12 sc_20 sc_21 sc_25 sc_30 sc_31 sc_36 sc_51 sc_52 sc_53
  batch_1 10293     0     0     0     0     0     0     0     0     0
  batch_2     0 12591  9658  7736     0     0     0     0     0     0
  batch_3     0     0     0     0 10195  8775  6006     0     0     0
  batch_4     0     0     0     0     0     0     0  6394  6241  8371

3. shr DA analysis

In [23]:
# subset only samples you want to compare
integrated.de <- subset(rc.integrated, subset = geno_trt %in% c("WT", "shr"))
In [24]:
DefaultAssay(integrated.de) <- "RNA"
merged <- as.SingleCellExperiment(integrated.de)
In [25]:
merged
class: SingleCellExperiment 
dim: 28429 72518 
metadata(0):
assays(2): counts logcounts
rownames(28429): AT1G01010 AT1G01020 ... AT4G07645 AT4G29600
rowData names(0):
colnames(72518): AAACCCAAGCCATTCA_1 AAACCCAAGCTGTTCA_1 ...
  TTTGTTGTCCTTTGAT_10 TTTGTTGTCTCGCTTG_10
colData names(138): orig.ident nCount_RNA ... sample.cell.p ident
reducedDimNames(5): PCA UMAP UMAP_50 UMAP_3D UMAP_2D
spikeNames(0):
altExpNames(0):
In [26]:
abundances <- table(merged$celltype.anno, merged$orig.ident) 
abundances <- unclass(abundances) 
abundances
A matrix: 12 × 8 of type int
sc_12sc_20sc_21sc_30sc_31sc_51sc_52sc_53
Quiescent Center 48 50 7 12 71 55 32 41
Stem Cell Niche 65 91 42 49 80 70 101 167
Columella13091116 608 866 829 77113221566
Lateral Root Cap13281446 93510821046113313831866
Atrichoblast13511556112015051052 705 7621302
Trichoblast1313186514571222 902 633 619 780
Cortex11461404108111481066 66310011264
Endodermis 664 947 9161361 959 666 95 95
Pericycle115911531100 895 777 524 135 183
Phloem 410 818 589 629 615 336 46 86
Xylem 404 821 668 548 507 288 115 188
Procambium109613241135 878 871 550 630 833
In [27]:
# Attaching some column metadata.
extra.info <- colData(merged)[match(colnames(abundances), merged$orig.ident),]
y.ab <- DGEList(abundances, samples=extra.info)
y.ab
$counts
A matrix: 12 × 8 of type int
sc_12sc_20sc_21sc_30sc_31sc_51sc_52sc_53
Quiescent Center 48 50 7 12 71 55 32 41
Stem Cell Niche 65 91 42 49 80 70 101 167
Columella13091116 608 866 829 77113221566
Lateral Root Cap13281446 93510821046113313831866
Atrichoblast13511556112015051052 705 7621302
Trichoblast1313186514571222 902 633 619 780
Cortex11461404108111481066 66310011264
Endodermis 664 947 9161361 959 666 95 95
Pericycle115911531100 895 777 524 135 183
Phloem 410 818 589 629 615 336 46 86
Xylem 404 821 668 548 507 288 115 188
Procambium109613241135 878 871 550 630 833
$samples
A data.frame: 8 × 141
grouplib.sizenorm.factorsorig.identnCount_RNAnFeature_RNAnCount_spliced_RNAnFeature_spliced_RNAnCount_unspliced_RNAnFeature_unspliced_RNAconsensus.time.grouptreatmenttimepointagerepdategenogeno_trtsample.cell.pident
<fct><dbl><dbl><chr><dbl><int><dbl><int><dbl><int><chr><chr><dbl><chr><int><chr><chr><fct><chr><fct>
sc_121102931sc_12117274244 0 0 0 0T3UntreatedNA5_dayNA2019-12-20WT WT sc_12,Atrichoblast 53
sc_201125911sc_20 43452318 0 0 0 0T7UntreatedNA5_dayNANA WT WT sc_20,Phloem 19
sc_211 96581sc_21406106449 0 0 0 0T4UntreatedNA5_dayNANA WT WT sc_21,Endodermis 38
sc_301101951sc_30 3770144933551325415278T8UntreatedNA5_dayNA2020-02-10WT WT sc_30,Cortex 12
sc_311 87751sc_31 8349295975282749821515T7UntreatedNA5_dayNA2020-02-10WT WT sc_31,Procambium 8
sc_511 63941sc_51129043395 0 0 0 0T9UntreatedNA5_dayNA2020-02-24WT WT sc_51,Xylem 21
sc_521 62411sc_52107103506 0 0 0 0T5UntreatedNA5_dayNA2020-02-24shr-2shrsc_52,Trichoblast 24
sc_531 83711sc_53 80983121 0 0 0 0T1UntreatedNA5_dayNA2020-02-24shr-2shrsc_53,Lateral Root Cap29
In [28]:
keep <- filterByExpr(y.ab, group=y.ab$samples$orig.ident)
y.ab <- y.ab[keep,]
summary(keep)
   Mode    TRUE 
logical      12 
In [29]:
(design <- model.matrix(~batch + factor(geno_trt), y.ab$samples))
A matrix: 8 × 5 of type dbl
(Intercept)batchbatch_2batchbatch_3batchbatch_4factor(geno_trt)shr
sc_1210000
sc_2011000
sc_2111000
sc_3010100
sc_3110100
sc_5110010
sc_5210011
sc_5310011
In [30]:
y.ab <- estimateDisp(y.ab, design, trend="none")
summary(y.ab$common.dispersion)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.03795 0.03795 0.03795 0.03795 0.03795 0.03795 
In [31]:
plotBCV(y.ab, cex=1)
In [32]:
fit.ab <- glmQLFit(y.ab, design, robust=TRUE, abundance.trend=FALSE)
summary(fit.ab$var.prior)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3863  0.3863  0.3863  0.3863  0.3863  0.3863 
In [33]:
summary(fit.ab$df.prior)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.5547  3.9487  3.9487  3.5185  3.9487  3.9487 
In [34]:
plotQLDisp(fit.ab, cex=1)
In [35]:
res <- glmQLFTest(fit.ab, coef=ncol(design))
summary(decideTests(res))
       factor(geno_trt)shr
Down                     4
NotSig                   6
Up                       2
In [36]:
topTags(res, n = nrow(res$table))
$table
A data.frame: 12 × 5
logFClogCPMFPValueFDR
<dbl><dbl><dbl><dbl><dbl>
Endodermis-2.9762783516.21968137.350323357.883747e-066.462717e-05
Phloem-2.5625209515.50018125.023440651.077120e-056.462717e-05
Pericycle-1.9131821616.24827 99.012964792.321556e-059.286225e-05
Xylem-1.1360199415.51238 34.442873816.357430e-041.907229e-03
Cortex 0.5862145616.90957 11.293837021.220667e-022.929601e-02
Columella 0.7257110316.89886 8.426388362.308058e-024.616116e-02
Lateral Root Cap 0.3268510317.17964 2.883555771.336122e-012.290495e-01
Atrichoblast 0.3328909016.97099 2.124913261.885834e-012.517830e-01
Stem Cell Niche 0.7256671613.27439 2.287335271.888372e-012.517830e-01
Procambium 0.2204835016.61436 1.381750342.785031e-013.342038e-01
Quiescent Center-0.7794572312.21389 0.314178596.085494e-016.638720e-01
Trichoblast-0.0419018516.84635 0.052749588.249601e-018.249601e-01
$adjust.method
'BH'
$comparison
'factor(geno_trt)shr'
$test
'glm'
In [37]:
shr_DA_results <- topTags(res, n = nrow(res$table))$table
In [38]:
shr_DA_results
A data.frame: 12 × 5
logFClogCPMFPValueFDR
<dbl><dbl><dbl><dbl><dbl>
Endodermis-2.9762783516.21968137.350323357.883747e-066.462717e-05
Phloem-2.5625209515.50018125.023440651.077120e-056.462717e-05
Pericycle-1.9131821616.24827 99.012964792.321556e-059.286225e-05
Xylem-1.1360199415.51238 34.442873816.357430e-041.907229e-03
Cortex 0.5862145616.90957 11.293837021.220667e-022.929601e-02
Columella 0.7257110316.89886 8.426388362.308058e-024.616116e-02
Lateral Root Cap 0.3268510317.17964 2.883555771.336122e-012.290495e-01
Atrichoblast 0.3328909016.97099 2.124913261.885834e-012.517830e-01
Stem Cell Niche 0.7256671613.27439 2.287335271.888372e-012.517830e-01
Procambium 0.2204835016.61436 1.381750342.785031e-013.342038e-01
Quiescent Center-0.7794572312.21389 0.314178596.085494e-016.638720e-01
Trichoblast-0.0419018516.84635 0.052749588.249601e-018.249601e-01

4. scr DA analysis

In [39]:
# subset only samples you want to compare
integrated.de <- subset(rc.integrated, subset = geno_trt %in% c("WT", "scr"))

DefaultAssay(integrated.de) <- "RNA"
merged <- as.SingleCellExperiment(integrated.de)

merged
class: SingleCellExperiment 
dim: 28429 71648 
metadata(0):
assays(2): counts logcounts
rownames(28429): AT1G01010 AT1G01020 ... AT4G07645 AT4G29600
rowData names(0):
colnames(71648): AAACCCAAGCCATTCA_1 AAACCCAAGCTGTTCA_1 ...
  TTTGTTGGTGGCTTGC_8 TTTGTTGTCTAAGCCA_8
colData names(138): orig.ident nCount_RNA ... sample.cell.p ident
reducedDimNames(5): PCA UMAP UMAP_50 UMAP_3D UMAP_2D
spikeNames(0):
altExpNames(0):
In [40]:
abundances <- table(merged$celltype.anno, merged$orig.ident) 
abundances <- unclass(abundances) 
abundances

# Attaching some column metadata.
extra.info <- colData(merged)[match(colnames(abundances), merged$orig.ident),]
y.ab <- DGEList(abundances, samples=extra.info)
y.ab

keep <- filterByExpr(y.ab, group=y.ab$samples$orig.ident)
y.ab <- y.ab[keep,]
summary(keep)
A matrix: 12 × 8 of type int
sc_12sc_20sc_21sc_25sc_30sc_31sc_36sc_51
Quiescent Center 48 50 7 45 12 71 36 55
Stem Cell Niche 65 91 42 56 49 80 110 70
Columella13091116 6081227 866 829 934 771
Lateral Root Cap13281446 93521341082104615031133
Atrichoblast135115561120115615051052 853 705
Trichoblast131318651457 8901222 9021061 633
Cortex114614041081 67311481066 273 663
Endodermis 664 947 916 2111361 959 150 666
Pericycle115911531100 247 895 777 206 524
Phloem 410 818 589 313 629 615 253 336
Xylem 404 821 668 282 548 507 143 288
Procambium109613241135 502 878 871 484 550
$counts
A matrix: 12 × 8 of type int
sc_12sc_20sc_21sc_25sc_30sc_31sc_36sc_51
Quiescent Center 48 50 7 45 12 71 36 55
Stem Cell Niche 65 91 42 56 49 80 110 70
Columella13091116 6081227 866 829 934 771
Lateral Root Cap13281446 93521341082104615031133
Atrichoblast135115561120115615051052 853 705
Trichoblast131318651457 8901222 9021061 633
Cortex114614041081 67311481066 273 663
Endodermis 664 947 916 2111361 959 150 666
Pericycle115911531100 247 895 777 206 524
Phloem 410 818 589 313 629 615 253 336
Xylem 404 821 668 282 548 507 143 288
Procambium109613241135 502 878 871 484 550
$samples
A data.frame: 8 × 141
grouplib.sizenorm.factorsorig.identnCount_RNAnFeature_RNAnCount_spliced_RNAnFeature_spliced_RNAnCount_unspliced_RNAnFeature_unspliced_RNAconsensus.time.grouptreatmenttimepointagerepdategenogeno_trtsample.cell.pident
<fct><dbl><dbl><chr><dbl><int><dbl><int><dbl><int><chr><chr><dbl><chr><int><chr><chr><fct><chr><fct>
sc_121102931sc_12117274244 0 0 0 0T3UntreatedNA5_dayNA2019-12-20WT WT sc_12,Atrichoblast 53
sc_201125911sc_20 43452318 0 0 0 0T7UntreatedNA5_dayNANA WT WT sc_20,Phloem 19
sc_211 96581sc_21406106449 0 0 0 0T4UntreatedNA5_dayNANA WT WT sc_21,Endodermis 38
sc_251 77361sc_25 94523249 0 0 0 0T5UntreatedNA5_dayNANA scr-4scrsc_25,Lateral Root Cap31
sc_301101951sc_30 3770144933551325415278T8UntreatedNA5_dayNA2020-02-10WT WT sc_30,Cortex 12
sc_311 87751sc_31 8349295975282749821515T7UntreatedNA5_dayNA2020-02-10WT WT sc_31,Procambium 8
sc_361 60061sc_36 57242386 0 0 0 0T7UntreatedNA5_dayNA2020-02-10scr-4scrsc_36,Endodermis 15
sc_511 63941sc_51129043395 0 0 0 0T9UntreatedNA5_dayNA2020-02-24WT WT sc_51,Xylem 21
   Mode    TRUE 
logical      12 
In [41]:
(design <- model.matrix(~batch + factor(geno_trt), y.ab$samples))
A matrix: 8 × 5 of type dbl
(Intercept)batchbatch_2batchbatch_3batchbatch_4factor(geno_trt)scr
sc_1210000
sc_2011000
sc_2111000
sc_2511001
sc_3010100
sc_3110100
sc_3610101
sc_5110010
In [42]:
y.ab <- estimateDisp(y.ab, design, trend="none")
summary(y.ab$common.dispersion)

plotBCV(y.ab, cex=1)

fit.ab <- glmQLFit(y.ab, design, robust=TRUE, abundance.trend=FALSE)
summary(fit.ab$var.prior)

summary(fit.ab$df.prior)

plotQLDisp(fit.ab, cex=1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.06526 0.06526 0.06526 0.06526 0.06526 0.06526 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.4215  0.4215  0.4215  0.4215  0.4215  0.4215 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.9332  3.6505  3.6505  3.3338  3.6505  3.6505 
In [43]:
res <- glmQLFTest(fit.ab, coef=ncol(design))
summary(decideTests(res))

topTags(res, n = nrow(res$table))

(scr_DA_results <- topTags(res, n = nrow(res$table))$table)
       factor(geno_trt)scr
Down                     5
NotSig                   5
Up                       2
$table
A data.frame: 12 × 5
logFClogCPMFPValueFDR
<dbl><dbl><dbl><dbl><dbl>
Endodermis-1.946839916.2783760.92401350.00013902600.0007259634
Pericycle-1.521518816.2987558.60748620.00015635360.0007259634
Lateral Root Cap 1.268225417.2758255.78276740.00018149080.0007259634
Xylem-1.048370415.5833333.01510900.00084231430.0025269429
Columella 0.931386316.7692921.96452410.00257077310.0061698554
Phloem-0.641658515.7294314.48193080.00733997700.0146799540
Cortex-0.829524316.6222110.06651440.01675264780.0287188247
Procambium-0.481334616.51235 5.46083440.05397026770.0809554015
Stem Cell Niche 0.898362513.07362 4.96460740.07089126080.0945216811
Atrichoblast 0.204461216.98999 1.34503450.28608499240.3433019908
Quiescent Center 0.833763812.28780 0.66719110.46061546780.5024896013
Trichoblast 0.174904616.98981 0.44696080.52629643560.5262964356
$adjust.method
'BH'
$comparison
'factor(geno_trt)scr'
$test
'glm'
A data.frame: 12 × 5
logFClogCPMFPValueFDR
<dbl><dbl><dbl><dbl><dbl>
Endodermis-1.946839916.2783760.92401350.00013902600.0007259634
Pericycle-1.521518816.2987558.60748620.00015635360.0007259634
Lateral Root Cap 1.268225417.2758255.78276740.00018149080.0007259634
Xylem-1.048370415.5833333.01510900.00084231430.0025269429
Columella 0.931386316.7692921.96452410.00257077310.0061698554
Phloem-0.641658515.7294314.48193080.00733997700.0146799540
Cortex-0.829524316.6222110.06651440.01675264780.0287188247
Procambium-0.481334616.51235 5.46083440.05397026770.0809554015
Stem Cell Niche 0.898362513.07362 4.96460740.07089126080.0945216811
Atrichoblast 0.204461216.98999 1.34503450.28608499240.3433019908
Quiescent Center 0.833763812.28780 0.66719110.46061546780.5024896013
Trichoblast 0.174904616.98981 0.44696080.52629643560.5262964356
In [44]:
y.ab2 <- calcNormFactors(y.ab)
y.ab2$samples$norm.factors

y.ab2 <- estimateDisp(y.ab2, design, trend="none")
fit.ab2 <- glmQLFit(y.ab2, design, robust=TRUE, abundance.trend=FALSE)
res2 <- glmQLFTest(fit.ab2, coef=ncol(design))
summary(decideTests(res2))
topTags(res, n = nrow(res$table))$table

(scr_DA_results_adj <- topTags(res, n = nrow(res$table))$table)
  1. 0.938282119936033
  2. 1.07092563921007
  3. 1.0763700673757
  4. 0.900568271322883
  5. 1.02034788507471
  6. 1.12562678551769
  7. 0.88093922949465
  8. 1.01470628783511
       factor(geno_trt)scr
Down                     3
NotSig                   7
Up                       2
A data.frame: 12 × 5
logFClogCPMFPValueFDR
<dbl><dbl><dbl><dbl><dbl>
Endodermis-1.946839916.2783760.92401350.00013902600.0007259634
Pericycle-1.521518816.2987558.60748620.00015635360.0007259634
Lateral Root Cap 1.268225417.2758255.78276740.00018149080.0007259634
Xylem-1.048370415.5833333.01510900.00084231430.0025269429
Columella 0.931386316.7692921.96452410.00257077310.0061698554
Phloem-0.641658515.7294314.48193080.00733997700.0146799540
Cortex-0.829524316.6222110.06651440.01675264780.0287188247
Procambium-0.481334616.51235 5.46083440.05397026770.0809554015
Stem Cell Niche 0.898362513.07362 4.96460740.07089126080.0945216811
Atrichoblast 0.204461216.98999 1.34503450.28608499240.3433019908
Quiescent Center 0.833763812.28780 0.66719110.46061546780.5024896013
Trichoblast 0.174904616.98981 0.44696080.52629643560.5262964356
A data.frame: 12 × 5
logFClogCPMFPValueFDR
<dbl><dbl><dbl><dbl><dbl>
Endodermis-1.946839916.2783760.92401350.00013902600.0007259634
Pericycle-1.521518816.2987558.60748620.00015635360.0007259634
Lateral Root Cap 1.268225417.2758255.78276740.00018149080.0007259634
Xylem-1.048370415.5833333.01510900.00084231430.0025269429
Columella 0.931386316.7692921.96452410.00257077310.0061698554
Phloem-0.641658515.7294314.48193080.00733997700.0146799540
Cortex-0.829524316.6222110.06651440.01675264780.0287188247
Procambium-0.481334616.51235 5.46083440.05397026770.0809554015
Stem Cell Niche 0.898362513.07362 4.96460740.07089126080.0945216811
Atrichoblast 0.204461216.98999 1.34503450.28608499240.3433019908
Quiescent Center 0.833763812.28780 0.66719110.46061546780.5024896013
Trichoblast 0.174904616.98981 0.44696080.52629643560.5262964356

5. Combine shr and scr DA results to plot

In [45]:
(shr_DA_results <- rownames_to_column(shr_DA_results, var = "celltype"))

shr_DA_results$geno <- rep("shr")

(scr_DA_results <- rownames_to_column(scr_DA_results, var = "celltype"))

scr_DA_results$geno <- rep("scr")

(DA_results <- bind_rows(shr_DA_results, scr_DA_results))

(DA_results <- mutate(DA_results, star = case_when(FDR <=0.001 ~ "***",  
                                          FDR <=0.01 ~ "**", 
                                          FDR <=0.05 ~ "*", 
                                          TRUE ~ " ")))
A data.frame: 12 × 6
celltypelogFClogCPMFPValueFDR
<chr><dbl><dbl><dbl><dbl><dbl>
Endodermis -2.9762783516.21968137.350323357.883747e-066.462717e-05
Phloem -2.5625209515.50018125.023440651.077120e-056.462717e-05
Pericycle -1.9131821616.24827 99.012964792.321556e-059.286225e-05
Xylem -1.1360199415.51238 34.442873816.357430e-041.907229e-03
Cortex 0.5862145616.90957 11.293837021.220667e-022.929601e-02
Columella 0.7257110316.89886 8.426388362.308058e-024.616116e-02
Lateral Root Cap 0.3268510317.17964 2.883555771.336122e-012.290495e-01
Atrichoblast 0.3328909016.97099 2.124913261.885834e-012.517830e-01
Stem Cell Niche 0.7256671613.27439 2.287335271.888372e-012.517830e-01
Procambium 0.2204835016.61436 1.381750342.785031e-013.342038e-01
Quiescent Center-0.7794572312.21389 0.314178596.085494e-016.638720e-01
Trichoblast -0.0419018516.84635 0.052749588.249601e-018.249601e-01
A data.frame: 12 × 6
celltypelogFClogCPMFPValueFDR
<chr><dbl><dbl><dbl><dbl><dbl>
Endodermis -1.946839916.2783760.92401350.00013902600.0007259634
Pericycle -1.521518816.2987558.60748620.00015635360.0007259634
Lateral Root Cap 1.268225417.2758255.78276740.00018149080.0007259634
Xylem -1.048370415.5833333.01510900.00084231430.0025269429
Columella 0.931386316.7692921.96452410.00257077310.0061698554
Phloem -0.641658515.7294314.48193080.00733997700.0146799540
Cortex -0.829524316.6222110.06651440.01675264780.0287188247
Procambium -0.481334616.51235 5.46083440.05397026770.0809554015
Stem Cell Niche 0.898362513.07362 4.96460740.07089126080.0945216811
Atrichoblast 0.204461216.98999 1.34503450.28608499240.3433019908
Quiescent Center 0.833763812.28780 0.66719110.46061546780.5024896013
Trichoblast 0.174904616.98981 0.44696080.52629643560.5262964356
A data.frame: 24 × 7
celltypelogFClogCPMFPValueFDRgeno
<chr><dbl><dbl><dbl><dbl><dbl><chr>
Endodermis -2.9762783516.21968137.350323357.883747e-066.462717e-05shr
Phloem -2.5625209515.50018125.023440651.077120e-056.462717e-05shr
Pericycle -1.9131821616.24827 99.012964792.321556e-059.286225e-05shr
Xylem -1.1360199415.51238 34.442873816.357430e-041.907229e-03shr
Cortex 0.5862145616.90957 11.293837021.220667e-022.929601e-02shr
Columella 0.7257110316.89886 8.426388362.308058e-024.616116e-02shr
Lateral Root Cap 0.3268510317.17964 2.883555771.336122e-012.290495e-01shr
Atrichoblast 0.3328909016.97099 2.124913261.885834e-012.517830e-01shr
Stem Cell Niche 0.7256671613.27439 2.287335271.888372e-012.517830e-01shr
Procambium 0.2204835016.61436 1.381750342.785031e-013.342038e-01shr
Quiescent Center-0.7794572312.21389 0.314178596.085494e-016.638720e-01shr
Trichoblast -0.0419018516.84635 0.052749588.249601e-018.249601e-01shr
Endodermis -1.9468398716.27837 60.924013531.390260e-047.259634e-04scr
Pericycle -1.5215188316.29875 58.607486161.563536e-047.259634e-04scr
Lateral Root Cap 1.2682254417.27582 55.782767361.814908e-047.259634e-04scr
Xylem -1.0483704215.58333 33.015108978.423143e-042.526943e-03scr
Columella 0.9313862816.76929 21.964524122.570773e-036.169855e-03scr
Phloem -0.6416584815.72943 14.481930837.339977e-031.467995e-02scr
Cortex -0.8295243316.62221 10.066514411.675265e-022.871882e-02scr
Procambium -0.4813346416.51235 5.460834425.397027e-028.095540e-02scr
Stem Cell Niche 0.8983624613.07362 4.964607417.089126e-029.452168e-02scr
Atrichoblast 0.2044611616.98999 1.345034472.860850e-013.433020e-01scr
Quiescent Center 0.8337638312.28780 0.667191094.606155e-015.024896e-01scr
Trichoblast 0.1749045716.98981 0.446960785.262964e-015.262964e-01scr
A data.frame: 24 × 8
celltypelogFClogCPMFPValueFDRgenostar
<chr><dbl><dbl><dbl><dbl><dbl><chr><chr>
Endodermis -2.9762783516.21968137.350323357.883747e-066.462717e-05shr***
Phloem -2.5625209515.50018125.023440651.077120e-056.462717e-05shr***
Pericycle -1.9131821616.24827 99.012964792.321556e-059.286225e-05shr***
Xylem -1.1360199415.51238 34.442873816.357430e-041.907229e-03shr**
Cortex 0.5862145616.90957 11.293837021.220667e-022.929601e-02shr*
Columella 0.7257110316.89886 8.426388362.308058e-024.616116e-02shr*
Lateral Root Cap 0.3268510317.17964 2.883555771.336122e-012.290495e-01shr
Atrichoblast 0.3328909016.97099 2.124913261.885834e-012.517830e-01shr
Stem Cell Niche 0.7256671613.27439 2.287335271.888372e-012.517830e-01shr
Procambium 0.2204835016.61436 1.381750342.785031e-013.342038e-01shr
Quiescent Center-0.7794572312.21389 0.314178596.085494e-016.638720e-01shr
Trichoblast -0.0419018516.84635 0.052749588.249601e-018.249601e-01shr
Endodermis -1.9468398716.27837 60.924013531.390260e-047.259634e-04scr***
Pericycle -1.5215188316.29875 58.607486161.563536e-047.259634e-04scr***
Lateral Root Cap 1.2682254417.27582 55.782767361.814908e-047.259634e-04scr***
Xylem -1.0483704215.58333 33.015108978.423143e-042.526943e-03scr**
Columella 0.9313862816.76929 21.964524122.570773e-036.169855e-03scr**
Phloem -0.6416584815.72943 14.481930837.339977e-031.467995e-02scr*
Cortex -0.8295243316.62221 10.066514411.675265e-022.871882e-02scr*
Procambium -0.4813346416.51235 5.460834425.397027e-028.095540e-02scr
Stem Cell Niche 0.8983624613.07362 4.964607417.089126e-029.452168e-02scr
Atrichoblast 0.2044611616.98999 1.345034472.860850e-013.433020e-01scr
Quiescent Center 0.8337638312.28780 0.667191094.606155e-015.024896e-01scr
Trichoblast 0.1749045716.98981 0.446960785.262964e-015.262964e-01scr
In [46]:
DA_results
A data.frame: 24 × 8
celltypelogFClogCPMFPValueFDRgenostar
<chr><dbl><dbl><dbl><dbl><dbl><chr><chr>
Endodermis -2.9762783516.21968137.350323357.883747e-066.462717e-05shr***
Phloem -2.5625209515.50018125.023440651.077120e-056.462717e-05shr***
Pericycle -1.9131821616.24827 99.012964792.321556e-059.286225e-05shr***
Xylem -1.1360199415.51238 34.442873816.357430e-041.907229e-03shr**
Cortex 0.5862145616.90957 11.293837021.220667e-022.929601e-02shr*
Columella 0.7257110316.89886 8.426388362.308058e-024.616116e-02shr*
Lateral Root Cap 0.3268510317.17964 2.883555771.336122e-012.290495e-01shr
Atrichoblast 0.3328909016.97099 2.124913261.885834e-012.517830e-01shr
Stem Cell Niche 0.7256671613.27439 2.287335271.888372e-012.517830e-01shr
Procambium 0.2204835016.61436 1.381750342.785031e-013.342038e-01shr
Quiescent Center-0.7794572312.21389 0.314178596.085494e-016.638720e-01shr
Trichoblast -0.0419018516.84635 0.052749588.249601e-018.249601e-01shr
Endodermis -1.9468398716.27837 60.924013531.390260e-047.259634e-04scr***
Pericycle -1.5215188316.29875 58.607486161.563536e-047.259634e-04scr***
Lateral Root Cap 1.2682254417.27582 55.782767361.814908e-047.259634e-04scr***
Xylem -1.0483704215.58333 33.015108978.423143e-042.526943e-03scr**
Columella 0.9313862816.76929 21.964524122.570773e-036.169855e-03scr**
Phloem -0.6416584815.72943 14.481930837.339977e-031.467995e-02scr*
Cortex -0.8295243316.62221 10.066514411.675265e-022.871882e-02scr*
Procambium -0.4813346416.51235 5.460834425.397027e-028.095540e-02scr
Stem Cell Niche 0.8983624613.07362 4.964607417.089126e-029.452168e-02scr
Atrichoblast 0.2044611616.98999 1.345034472.860850e-013.433020e-01scr
Quiescent Center 0.8337638312.28780 0.667191094.606155e-015.024896e-01scr
Trichoblast 0.1749045716.98981 0.446960785.262964e-015.262964e-01scr
In [47]:
write.csv(DA_results, "./supp_data/shr_scr_DA_celltype.csv", row.names=F)
In [48]:
FC_wide <- select(DA_results, geno, , celltype, logFC) %>%
  pivot_wider(names_from = geno, values_from = logFC) %>% mutate(WT=rep(0))

FC_wide_m <- as.matrix(FC_wide[, 2:ncol(FC_wide)])
rownames(FC_wide_m) <- FC_wide$celltype

FC_wide_m

star.wide <- select(DA_results, geno, , celltype, star) %>%
  pivot_wider(names_from = geno, values_from = star) %>%
mutate(WT=rep(" "))

star.wide[is.na(star.wide)] <- " "

star.wide_m <- as.matrix(star.wide[, 2:ncol(star.wide)])
rownames(star.wide_m) <- star.wide$celltype
FC_wide_m
star.wide_m
A matrix: 12 × 3 of type dbl
shrscrWT
Endodermis-2.97627835-1.94683990
Phloem-2.56252095-0.64165850
Pericycle-1.91318216-1.52151880
Xylem-1.13601994-1.04837040
Cortex 0.58621456-0.82952430
Columella 0.72571103 0.93138630
Lateral Root Cap 0.32685103 1.26822540
Atrichoblast 0.33289090 0.20446120
Stem Cell Niche 0.72566716 0.89836250
Procambium 0.22048350-0.48133460
Quiescent Center-0.77945723 0.83376380
Trichoblast-0.04190185 0.17490460
A matrix: 12 × 3 of type dbl
shrscrWT
Endodermis-2.97627835-1.94683990
Phloem-2.56252095-0.64165850
Pericycle-1.91318216-1.52151880
Xylem-1.13601994-1.04837040
Cortex 0.58621456-0.82952430
Columella 0.72571103 0.93138630
Lateral Root Cap 0.32685103 1.26822540
Atrichoblast 0.33289090 0.20446120
Stem Cell Niche 0.72566716 0.89836250
Procambium 0.22048350-0.48133460
Quiescent Center-0.77945723 0.83376380
Trichoblast-0.04190185 0.17490460
A matrix: 12 × 3 of type chr
shrscrWT
Endodermis******
Phloem****
Pericycle******
Xylem** **
Cortex* *
Columella* **
Lateral Root Cap ***
Atrichoblast
Stem Cell Niche
Procambium
Quiescent Center
Trichoblast
In [49]:
options(repr.plot.width=6, repr.plot.height=12)

vert_Cell_hm <- Heatmap(FC_wide_m, 
                name = "Mutant/WT Cell type", 
                col = colorRamp2(c(-4, 0, 4), c("blue", "beige", "#e31a1c")), 
                column_title = " ", 
                cluster_rows = T,
                cluster_columns = F, 
                clustering_distance_columns = "pearson", 
                clustering_distance_rows = "pearson",
                use_raster= FALSE, 
                show_column_names = TRUE, 
                show_row_names = TRUE, 
                show_row_dend = TRUE, 
                show_column_dend = TRUE, 
                        heatmap_legend_param = list(title_position="topcenter", 
                                               color_bar = "continuous", 
                                               legend_direction = "horizontal"), 
                cell_fun = function(j, i, x, y, width, height, fill) 
                  {grid.text(sprintf("%s", star.wide_m[i, j]), x, y, gp = gpar(fontsize = 10))}) 

draw(vert_Cell_hm, padding = unit(c(10, 10, 10, 10), "mm"), heatmap_legend_side = "top")
In [50]:
options(repr.plot.width=3.6, repr.plot.height=9)
	



FC_wide_m_2 <- FC_wide_m[c("Endodermis", 
                           "Cortex", 
                           "Pericycle", 
                           "Phloem", 
                           "Xylem",  
                           "Procambium", 
                           "Quiescent Center", 
                           "Trichoblast", 
                           "Atrichoblast", 
                           "Stem Cell Niche", 
                           "Columella", 
                           "Lateral Root Cap") , 
                         c("shr", "scr")]

star.wide_m_2 <- star.wide_m[c("Endodermis", 
                           "Cortex", 
                           "Pericycle", 
                           "Phloem", 
                           "Xylem",  
                           "Procambium", 
                               "Quiescent Center", 
                           "Trichoblast", 
                           "Atrichoblast", 
                            "Stem Cell Niche", 
                           "Columella", 
                           "Lateral Root Cap") , 
                         c("shr", "scr")]

col_lab_it <- c( 
              expression(italic("shr ")),
                expression(italic("scr ")))

vert_Cell_hm2 <- Heatmap(FC_wide_m_2, 
                name = "Cell Type logFC \n (mutant/WT)", 
                col = colorRamp2(c(-3, 0, 3), c("blue", "beige", "#e31a1c")), 
                column_title = " ", 
                cluster_rows = F,
                cluster_columns = F, 
                clustering_distance_columns = "pearson", 
                clustering_distance_rows = "pearson",
                use_raster= FALSE, 
                show_column_names = TRUE, 
                show_row_names = TRUE, 
                show_row_dend = TRUE, 
                show_column_dend = TRUE, 
                         column_labels=col_lab_it,
                         column_names_gp = gpar(fontsize = 16),
                         row_names_gp = gpar(fontsize = 16),
                        heatmap_legend_param = list(title_position="topcenter", 
                                               color_bar = "continuous", 
                                               legend_direction = "horizontal", 
                                                   labels_gp = gpar(fontsize = 16), 
                                                   title_gp = gpar(fontsize = 16, font = 2), 
                                                   legend_height = unit(4, "cm")), 
                cell_fun = function(j, i, x, y, width, height, fill) 
                  {grid.text(sprintf("%s", star.wide_m_2[i, j]), x, y, gp = gpar(fontsize = 14))}) 

draw(vert_Cell_hm2, padding = unit(c(10, 10, 10, 10), "mm"), heatmap_legend_side = "top")


hm_final <- grid.grabExpr(draw(vert_Cell_hm2, padding = unit(c(10, 10, 10, 10), "mm"), heatmap_legend_side = "top"))
In [51]:
FC_wide_m_2 <- t(FC_wide_m_2)
rownames(FC_wide_m_2) <- c("shr-2", "scr-4")
star.wide_m_2 <- t(star.wide_m_2)
rownames(star.wide_m_2) <- c("shr-2", "scr-4")
In [52]:
FC_wide_m_2
star.wide_m_2
A matrix: 2 × 12 of type dbl
EndodermisCortexPericyclePhloemXylemProcambiumQuiescent CenterTrichoblastAtrichoblastStem Cell NicheColumellaLateral Root Cap
shr-2-2.976278 0.5862146-1.913182-2.5625210-1.13602 0.2204835-0.7794572-0.041901850.33289090.72566720.72571100.326851
scr-4-1.946840-0.8295243-1.521519-0.6416585-1.04837-0.4813346 0.8337638 0.174904570.20446120.89836250.93138631.268225
A matrix: 2 × 12 of type chr
EndodermisCortexPericyclePhloemXylemProcambiumQuiescent CenterTrichoblastAtrichoblastStem Cell NicheColumellaLateral Root Cap
shr-2************ *
scr-4******** ** *****
In [53]:
options(repr.plot.width=8, repr.plot.height=4)

hor_Cell_hm2 <- Heatmap(FC_wide_m_2, 
                name = "Cell Type logFC \n (mutant/WT)", 
                col = colorRamp2(c(-3, 0, 3), c("blue", "beige", "#e31a1c")), 
                column_title = " ", 
                cluster_rows = F,
                cluster_columns = F, 
                clustering_distance_columns = "pearson", 
                clustering_distance_rows = "pearson",
                use_raster= FALSE, 
                show_column_names = TRUE, 
                show_row_names = TRUE, 
                show_row_dend = TRUE, 
                        row_names_side = "left",
                show_column_dend = TRUE, 
                        column_names_rot = 45,
                        heatmap_legend_param = list(title_position="topcenter", 
                                               color_bar = "continuous", 
                                               legend_direction = "horizontal"), 
                        row_names_gp = gpar(fontsize = 14, font = 3),
                cell_fun = function(j, i, x, y, width, height, fill) 
                  {grid.text(sprintf("%s", star.wide_m_2[i, j]), x, y, gp = gpar(fontsize = 14))}) 

draw(hor_Cell_hm2, padding = unit(c(10, 20, 10, 10), "mm"), heatmap_legend_side = "top")

pdf("./supp_data/shr_scr_horizontal_DA.pdf", width = 8, height = 4)
draw(hor_Cell_hm2, padding = unit(c(10, 20, 10, 10), "mm"), heatmap_legend_side = "top")
dev.off()
png: 2

6. Label transfer classification/prediction score

In [54]:
anno_scores <- data.frame(geno=rc.integrated$geno_trt, 
           prediction.score.Endodermis=rc.integrated$prediction.score.Endodermis, 
           prediction.score.Cortex=rc.integrated$prediction.score.Cortex, 
           annotation=rc.integrated$celltype.anno, 
                         timezone=rc.integrated$time.anno, 
                         consensus.time=rc.integrated$consensus.time.group)
In [55]:
anno_scores %>% filter(annotation %in% c("Cortex", "Endodermis")) %>% filter(is.na(prediction.score.Cortex)==F)
A data.frame: 15783 × 6
genoprediction.score.Endodermisprediction.score.Cortexannotationtimezoneconsensus.time
<fct><dbl><dbl><fct><fct><fct>
WT0.000000001.00000000Cortex ElongationT6
WT0.000000001.00000000Cortex ElongationT6
WT1.000000000.00000000EndodermisElongationT4
WT1.000000000.00000000EndodermisMaturationT6
WT0.017169840.98283016Cortex MaturationT7
WT0.000000001.00000000Cortex MaturationT9
WT1.000000000.00000000EndodermisMeristem T1
WT0.000000001.00000000Cortex ElongationT9
WT0.832070820.06699925EndodermisMeristem T2
WT1.000000000.00000000EndodermisMeristem T2
WT0.000000001.00000000Cortex ElongationT7
WT1.000000000.00000000EndodermisElongationT5
WT0.000000001.00000000Cortex MaturationT9
WT0.000000001.00000000Cortex ElongationT5
WT0.796958220.00000000EndodermisMaturationT5
WT0.000000001.00000000Cortex MaturationT7
WT0.721752030.00000000EndodermisMaturationT5
WT0.000000001.00000000Cortex ElongationT7
WT0.000000001.00000000Cortex ElongationT7
WT1.000000000.00000000EndodermisMaturationT6
WT1.000000000.00000000EndodermisElongationT4
WT0.000000001.00000000Cortex MaturationT8
WT0.000000001.00000000Cortex ElongationT5
WT0.000000001.00000000Cortex MaturationT9
WT1.000000000.00000000EndodermisElongationT2
WT1.000000000.00000000EndodermisMeristem T2
WT1.000000000.00000000EndodermisMeristem T2
WT0.000000001.00000000Cortex ElongationT8
WT1.000000000.00000000EndodermisElongationT9
WT0.000000001.00000000Cortex ElongationT4
shr0.0000000001.0000000Cortex Meristem T2
shr0.0051773250.9948227Cortex MaturationT8
shr0.0051097070.8627821Cortex Meristem T1
shr0.0000000001.0000000Cortex ElongationT5
shr0.0846397950.9153602Cortex ElongationT7
shr0.0000000001.0000000Cortex ElongationT6
shr0.6609778120.2107672EndodermisMeristem T3
shr0.0587781770.9412218Cortex MaturationT7
shr0.0930536800.9069463Cortex ElongationT7
shr0.0186362030.9813638Cortex ElongationT4
shr0.0440875300.9529055Cortex ElongationT8
shr0.0000000001.0000000Cortex MaturationT8
shr0.0000000001.0000000Cortex MaturationT8
shr0.0000000001.0000000Cortex ElongationT6
shr0.0209041040.9763210Cortex ElongationT8
shr0.0000000001.0000000Cortex MaturationT8
shr0.0083776650.9194263Cortex ElongationT3
shr0.6913945610.3043475EndodermisElongationT5
shr0.0000000001.0000000Cortex ElongationT6
shr0.0147602190.9852398Cortex MaturationT8
shr0.1835252460.6478574Cortex ElongationT8
shr0.0000000001.0000000Cortex ElongationT3
shr0.0309155100.9690845Cortex ElongationT8
shr0.2374442080.7625558Cortex ElongationT7
shr0.0011132330.9988868Cortex Meristem T2
shr0.0000000000.9734162Cortex Meristem T2
shr0.0070150520.9787497Cortex MaturationT9
shr0.0000000001.0000000Cortex MaturationT8
shr0.0194022990.9805977Cortex ElongationT8
shr0.0000000001.0000000Cortex MaturationT7
In [56]:
pred_plot <- anno_scores %>% filter(annotation %in% c("Cortex", "Endodermis")) %>% filter(is.na(prediction.score.Cortex)==F) %>%
ggplot(aes(x=prediction.score.Endodermis, y=prediction.score.Cortex, color=annotation)) + geom_point(alpha=0.1)
In [57]:
library(ggridges)
In [58]:
options(repr.plot.width=8, repr.plot.height=6)
pred_plot <- anno_scores %>% filter(annotation %in% c("Cortex", "Endodermis")) %>% filter(is.na(prediction.score.Cortex)==F) %>% 
mutate(Endodermis_Cortex_diff=prediction.score.Cortex-prediction.score.Endodermis)  %>%
ggplot(aes(y=factor(geno, levels=c("scr", "shr", "WT")), x=Endodermis_Cortex_diff, fill=geno)) + geom_density_ridges(rel_min_height = 0.00) 

mylabels <- c( 
              expression(italic("scr-4")),
                expression(italic("shr-2")), "WT")

pred_plot <- pred_plot + 
theme_cowplot(font_size = 20) + 
xlab("Cortex Prediction score - Endodermis Prediction score") + 
ylab("") + 
scale_fill_viridis_d(labels = rev(mylabels)) + 
scale_x_continuous(limits=c(-1.2, 1.2)) + 
scale_y_discrete(labels = mylabels)
pred_plot
Picking joint bandwidth of 0.102

In [59]:
options(repr.plot.width=8, repr.plot.height=10)
pred_plot <- pred_plot + facet_wrap(~ factor(timezone, levels=c("Maturation","Elongation","Meristem")), ncol = 1) + theme(legend.title = element_blank())
pred_plot
Picking joint bandwidth of 0.13

Picking joint bandwidth of 0.0907

Picking joint bandwidth of 0.0809

In [60]:
# make triangles 
options(repr.plot.width=8, repr.plot.height=2)

d=data.frame(x=c(1,1,16, 16,16,1), y=c(1,2,1.5, 1.7,2.7,2.2), t=c('a', 'a', 'a', 'b', 'b', 'b'))
d
triangle_label <- ggplot() +
geom_polygon(data=d, mapping=aes(x=x, y=y, group=t, fill=t)) + annotate("text", x = 3.5, y = 1.5, label = "Endodermis", size = 8, color="white") + 
annotate("text", x = 14.5, y = 2.2, label = "Cortex", size = 8, color="white") + 
scale_fill_manual(values=(c("#0000FF", "#82B6FF"))) + theme_nothing() + scale_x_continuous(limits=c(0.1, 17))

triangle_label
A data.frame: 6 × 3
xyt
<dbl><dbl><fct>
11.0a
12.0a
161.5a
161.7b
162.7b
12.2b
In [61]:
options(repr.plot.width=9, repr.plot.height=11)

(pred_tri <- plot_grid(pred_plot, triangle_label, ncol=1, rel_heights = c(1, .2)))

ggsave("./supp_data/WT_shr_scr_prediction_ridgeplot.pdf", width=9, height=10)

ggsave("./supp_data/WT_shr_scr_prediction_ridgeplot_taller.pdf", width=9, height=11)

ggsave("./supp_data/WT_shr_scr_prediction_ridgeplot_tallest.pdf", width=9, height=12)
Picking joint bandwidth of 0.13

Picking joint bandwidth of 0.0907

Picking joint bandwidth of 0.0809

In [ ]:

In [62]:
options(repr.plot.width=15, repr.plot.height=6)
time_plt <- DimPlot(rc.integrated, 
        group.by = "time.anno", 
        order = c("Meristem","Elongation","Maturation"),
        cols = c("#1A237E", "#42B3D5", "#DCEDC8"), split.by="geno_trt")
time_plt
In [63]:
options(repr.plot.width=12, repr.plot.height=12)

cell_plt <- DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", split.by="geno_trt", cols=color)


(cell_time_plt <- plot_grid(cell_plt, time_plt, nrow=2))
In [64]:
options(repr.plot.width=20, repr.plot.height=12)

plot_grid(cell_time_plt, pred_tri, ncol=2, rel_widths = c(2,1.2))
In [ ]:

In [65]:
# swtich direction
In [66]:
options(repr.plot.width=8, repr.plot.height=6)
pred_plot <- anno_scores %>% filter(annotation %in% c("Cortex", "Endodermis")) %>% filter(is.na(prediction.score.Cortex)==F) %>% 
mutate(Cortex__Endo_diff=prediction.score.Endodermis-prediction.score.Cortex)  %>%
ggplot(aes(y=factor(geno, levels=c("shr", "scr", "WT")), x=Cortex__Endo_diff, fill=geno)) + geom_density_ridges(rel_min_height = 0.00) 

mylabels <- c( 
              expression(italic("shr")),
                expression(italic("scr")), "WT")

pred_plot <- pred_plot + 
theme_cowplot() + 
xlab("Endodermis Prediction score - Cortex Prediction score") + 
ylab("") + 
scale_fill_viridis_d(labels = rev(mylabels)) + 
scale_x_continuous(limits=c(-1.2, 1.2)) + 
scale_y_discrete(labels = mylabels) + theme(legend.position = "none")
pred_plot
Picking joint bandwidth of 0.102

In [ ]:

In [67]:
options(repr.plot.width=8, repr.plot.height=10)
pred_plot <- pred_plot + facet_wrap(~ factor(timezone, levels=c("Maturation","Elongation","Meristem")), ncol = 1) + theme(legend.title = element_blank())
pred_plot
Picking joint bandwidth of 0.13

Picking joint bandwidth of 0.0907

Picking joint bandwidth of 0.0809

In [68]:
# make triangles 
options(repr.plot.width=8, repr.plot.height=2)

d=data.frame(x=c(1,1,16, 16,16,1), y=c(1,2,1.5, 1.7,2.7,2.2), t=c('a', 'a', 'a', 'b', 'b', 'b'))
d
triangle_label <- ggplot() +
geom_polygon(data=d, mapping=aes(x=x, y=y, group=t, fill=t)) + annotate("text", x = 2.5, y = 1.5, label = "Cortex", size = 7, color="white") + 
annotate("text", x = 13.8, y = 2.2, label = "Endodermis", size = 7, color="white") + 
scale_fill_manual(values=(c("blue", "red"))) + theme_nothing() + scale_x_continuous(limits=c(0.1, 17))

triangle_label
A data.frame: 6 × 3
xyt
<dbl><dbl><fct>
11.0a
12.0a
161.5a
161.7b
162.7b
12.2b
In [69]:
options(repr.plot.width=8, repr.plot.height=10)

(pred_tri <- plot_grid(pred_plot, triangle_label, ncol=1, rel_heights = c(1, .15), align="v"))
Picking joint bandwidth of 0.13

Picking joint bandwidth of 0.0907

Picking joint bandwidth of 0.0809

In [70]:
rc.integrated$Endodermis_Cortex_diff <- rc.integrated$prediction.score.Cortex - rc.integrated$prediction.score.Endodermis
In [71]:
options(repr.plot.width=15, repr.plot.height=6)

FeaturePlot(rc.integrated, features = "Endodermis_Cortex_diff",
    cols = c("blue", "grey", "red"), label=F, repel=F, pt.size = 0.5, split.by="geno_trt")
In [72]:
options(repr.plot.width=15, repr.plot.height=18)
# Visualize co-expression of two features simultaneously
FeaturePlot(rc.integrated, features = c("prediction.score.Cortex", "prediction.score.Endodermis"), blend = TRUE, split.by="geno_trt", cols = c("grey","blue", "red"))
In [73]:
WT_ground <- readRDS("./supp_data/Ground_Tissue_Atlas.rds")
In [74]:
scr_ground <- readRDS("./supp_data/Ground_Tissue_scr.rds")
In [75]:
shr_ground <- readRDS("./supp_data/Ground_Tissue_shr.rds")
In [76]:
options(repr.plot.width=15, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
FeaturePlot(scr_ground, features = c("prediction.score.Cortex", "prediction.score.Endodermis"), blend = TRUE, cols = c("grey","blue", "red"), min.cutoff = 0, max.cutoff = 1)
In [77]:
summary(scr_ground$prediction.score.Cortex)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.02845 0.74804 0.56686 0.95337 1.00000 
In [78]:
options(repr.plot.width=15, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
FeaturePlot(shr_ground, features = c("prediction.score.Cortex", "prediction.score.Endodermis"), blend = TRUE, cols = c("grey","blue", "red"), min.cutoff = 0, max.cutoff = 1)
In [79]:
options(repr.plot.width=6, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
(scr_e <- FeaturePlot(scr_ground, features = c("prediction.score.Endodermis"), cols = c("grey", "red"), min.cutoff = 0, max.cutoff = 1) + ggtitle("scr \n Prediction Score \n Endodermis"))
In [80]:
options(repr.plot.width=6, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
(scr_c <- FeaturePlot(scr_ground, features = c("prediction.score.Cortex"), cols = c("grey", "blue"), min.cutoff = 0, max.cutoff = 1) + ggtitle("scr \n Prediction Score \n Cortex"))
In [81]:
scr_time_plt <- DimPlot(scr_ground, 
        group.by = "time.anno", 
        order = c("Meristem","Elongation","Maturation"),
        cols = c("#1A237E", "#42B3D5", "#DCEDC8")) + ggtitle("scr \n Developmental Stage") + theme(plot.title = element_text(hjust = 0.5))
scr_time_plt
In [82]:
options(repr.plot.width=6, repr.plot.height=8)
scr_ct_plt <-DimPlot(scr_ground, reduction = "umap", group.by = "consensus.time.group", cols=brewer.pal(10,"Spectral")) + ggtitle("scr \n Consensus Time") + theme(plot.title = element_text(hjust = 0.5))
scr_ct_plt
In [ ]:

In [83]:
options(repr.plot.width=6, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
(shr_e <- FeaturePlot(shr_ground, features = c("prediction.score.Endodermis"), cols = c("grey", "red"), min.cutoff = 0, max.cutoff = 1) + ggtitle("shr \n Prediction Score \n Endodermis"))
In [84]:
options(repr.plot.width=6, repr.plot.height=8)
# Visualize co-expression of two features simultaneously
(shr_c <- FeaturePlot(shr_ground, features = c("prediction.score.Cortex"), cols = c("grey", "blue"), min.cutoff = 0, max.cutoff = 1) + ggtitle("shr \n Prediction Score \n Cortex"))
In [85]:
shr_time_plt <- DimPlot(shr_ground, 
        group.by = "time.anno", 
        order = c("Meristem","Elongation","Maturation"),
        cols = c("#1A237E", "#42B3D5", "#DCEDC8")) + ggtitle("shr \n Developmental Stage") + theme(plot.title = element_text(hjust = 0.5))
shr_time_plt
In [86]:
options(repr.plot.width=6, repr.plot.height=8)
shr_ct_plt <-DimPlot(shr_ground, reduction = "umap", group.by = "consensus.time.group", cols=brewer.pal(10,"Spectral")) + ggtitle("shr \n Consensus Time") + theme(plot.title = element_text(hjust = 0.5))
shr_ct_plt
In [87]:
options(repr.plot.width=18, repr.plot.height=16)
(pred_dim_plots <- plot_grid(scr_time_plt, scr_c, scr_e, shr_time_plt, shr_c, shr_e, ncol=3, align="hv"))
In [88]:
options(repr.plot.width=24, repr.plot.height=12)

plot_grid(pred_dim_plots, pred_tri, ncol=2, rel_widths = c(2.5,1.2))

ggsave("./supp_data/scr_shr_pred_dim.pdf", width=30, height=12)
In [89]:
options(repr.plot.width=8, repr.plot.height=6)
pred_plot <- anno_scores %>% filter(annotation %in% c("Cortex", "Endodermis")) %>% filter(is.na(prediction.score.Cortex)==F) %>% 
mutate(Cortex__Endo_diff=prediction.score.Endodermis-prediction.score.Cortex)  %>%
ggplot(aes(y=factor(geno, levels=c("shr", "scr", "WT")), x=Cortex__Endo_diff, fill=geno)) + geom_density_ridges(rel_min_height = 0.00) 

mylabels <- c( 
              expression(italic("shr")),
                expression(italic("scr")), "WT")

pred_plot <- pred_plot + 
theme_cowplot() + 
xlab("Endodermis Prediction score - Cortex Prediction score") + 
ylab("") + 
scale_fill_viridis_d(labels = rev(mylabels)) + 
scale_x_continuous(limits=c(-1.2, 1.2)) + 
scale_y_discrete(labels = mylabels) + theme(legend.position = "none")
pred_plot
Picking joint bandwidth of 0.102

In [90]:
options(repr.plot.width=8, repr.plot.height=15)

con_time_order <- c("T0", "T1", "T2", "T3", "T4", "T5", "T6", "T7", "T8", "T9")
con_time_order <- rev(con_time_order)


pred_plot <- pred_plot + facet_wrap(~ factor(consensus.time, levels=con_time_order), ncol = 1) + theme(legend.title = element_blank())
pred_plot
Picking joint bandwidth of 0.146

Picking joint bandwidth of 0.14

Picking joint bandwidth of 0.0599

Picking joint bandwidth of 0.117

Picking joint bandwidth of 0.128

Picking joint bandwidth of 0.0892

Picking joint bandwidth of 0.0854

Picking joint bandwidth of 0.185

Picking joint bandwidth of 0.15

Picking joint bandwidth of 0.231

In [91]:
options(repr.plot.width=8, repr.plot.height=20)

(pred_tri <- plot_grid(pred_plot, triangle_label, ncol=1, rel_heights = c(1, .1)))
Picking joint bandwidth of 0.146

Picking joint bandwidth of 0.14

Picking joint bandwidth of 0.0599

Picking joint bandwidth of 0.117

Picking joint bandwidth of 0.128

Picking joint bandwidth of 0.0892

Picking joint bandwidth of 0.0854

Picking joint bandwidth of 0.185

Picking joint bandwidth of 0.15

Picking joint bandwidth of 0.231

In [ ]:

In [92]:
options(repr.plot.width=24, repr.plot.height=16)
(pred_dim_plots <- plot_grid(scr_time_plt, scr_ct_plt, scr_c, scr_e, shr_time_plt, shr_ct_plt, shr_c, shr_e, ncol=4, align="hv"))
In [93]:
options(repr.plot.width=30, repr.plot.height=15)

plot_grid(pred_dim_plots, pred_tri, ncol=2, rel_widths = c(3.5,1.2))

ggsave("./supp_data/scr_shr_pred_dim_with_consensus_time.pdf", width=30, height=15)
In [ ]: